Regresión Lineal Múltiple

TP Final





Alumnos

Nicolás Dominutti, Carlos Suarez Gurruchaga, Hernán Telechea







Primer acercamiento y correcciones

SegVial <- read_excel("SegVial.xlsx")

(head(SegVial,3))
## # A tibble: 3 × 20
##   `Codigo Pais` Pais   Ciudad   Poblacion DenPob ArCiclista ArBajaVel PMPeatones
##   <chr>         <chr>  <chr>        <dbl>  <dbl> <chr>      <chr>     <chr>     
## 1 ES            Spain  Barcelo…   1620343  16014 0.15       0.775     32.0      
## 2 ES            Spain  Madrid     3223334   5332 0.067      0.065     34.0      
## 3 FR            France Bordeaux    254436   5104 0.106      0.239     21.0      
## # … with 12 more variables: PMCiclistas <chr>, PMTPublico <chr>,
## #   PMVMotor <chr>, PeatAuto <dbl>, CicAuto <dbl>, V2RMSM <dbl>,
## #   V2RMAuto <dbl>, AutoSM <dbl>, AutoAuto <dbl>, Precipitacion <chr>,
## #   Temp <chr>, PBI <dbl>
glimpse(SegVial)
## Rows: 24
## Columns: 20
## $ `Codigo Pais` <chr> "ES", "ES", "FR", "FR", "FR", "FR", "FR", "FR", "FR", "F…
## $ Pais          <chr> "Spain", "Spain", "France", "France", "France", "France"…
## $ Ciudad        <chr> "Barcelona", "Madrid", "Bordeaux", "Lille", "Lyon", "Mar…
## $ Poblacion     <dbl> 1620343, 3223334, 254436, 232787, 516092, 863310, 285121…
## $ DenPob        <dbl> 16014, 5332, 5104, 6687, 10758, 3564, 4997, 4704, 4586, …
## $ ArCiclista    <chr> "0.15", "0.067", "0.106", "0.098", "0.109", "0.044", "0.…
## $ ArBajaVel     <chr> "0.775", "0.065", "0.239", "0.634", "0.284", "0.118", "0…
## $ PMPeatones    <chr> "32.0", "34.0", "21.0", "32.0", "34.0", "34.0", "26.0", …
## $ PMCiclistas   <chr> "2.0", "0.46", "3.0", "2.0", "2.0", "1.0", "2.0", "5.0",…
## $ PMTPublico    <chr> "39.0", "24.35", "9.0", "10.0", "19.0", "11.0", "8.0", "…
## $ PMVMotor      <chr> "27.0", "39.89", "67.0", "56.0", "45.0", "54.0", "64.0",…
## $ PeatAuto      <dbl> 22, 201, 13, 5, 31, 56, 21, 13, 20, 58, 5, 12, 153, 60, …
## $ CicAuto       <dbl> 5, 18, 5, 1, 11, 9, 4, 7, 2, 14, 8, 2, 54, 14, 28, 15, 2…
## $ V2RMSM        <dbl> 29, 95, 2, 2, 12, 45, 8, 6, 6, 48, 4, 5, 6, 8, 1, 1, 2, …
## $ V2RMAuto      <dbl> 96, 275, 14, 7, 39, 146, 21, 13, 50, 132, 8, 29, 55, 11,…
## $ AutoSM        <dbl> 0, 44, 3, 2, 8, 30, 9, 1, 5, 12, 4, 5, 20, 18, 6, 7, 11,…
## $ AutoAuto      <dbl> 4, 72, 1, 6, 21, 48, 12, 6, 6, 24, 3, 12, 122, 31, 11, 1…
## $ Precipitacion <chr> "565.0", "420.9", "944.1", "742.5", "831.9", "515.4", "6…
## $ Temp          <chr> "18.2", "15.0", "13.8", "10.8", "12.5", "15.5", "15.1", …
## $ PBI           <dbl> 31070.56, 35290.16, 35464.54, 31291.93, 48140.02, 36231.…

Tenemos 20 variables analizadas con 24 datos cada una, no se observan valores NaNs. 11 de esas variables aparecen como categóricas, el resto como cuantitativas. Sin embargo, al observar la tabla de datos, solo 3 deberían ser categóricas. Las demás aparecen como tal por tener decimales. Las vamos a transformar.

SegVial[,c(6, 7, 8, 9, 10, 11, 18, 19)] = lapply(SegVial[,c(6, 7, 8, 9, 10, 11, 18, 19)], function(x) as.numeric(x))
summary(SegVial)
##  Codigo Pais            Pais              Ciudad            Poblacion      
##  Length:24          Length:24          Length:24          Min.   : 232787  
##  Class :character   Class :character   Class :character   1st Qu.: 432558  
##  Mode  :character   Mode  :character   Mode  :character   Median : 542400  
##                                                           Mean   : 976996  
##                                                           3rd Qu.: 932826  
##                                                           Max.   :3600203  
##      DenPob        ArCiclista        ArBajaVel        PMPeatones   
##  Min.   : 1417   Min.   :0.02700   Min.   :0.0240   Min.   : 4.00  
##  1st Qu.: 3223   1st Qu.:0.06425   1st Qu.:0.1003   1st Qu.:11.50  
##  Median : 4212   Median :0.10200   Median :0.2380   Median :26.50  
##  Mean   : 5483   Mean   :0.11675   Mean   :0.2974   Mean   :23.98  
##  3rd Qu.: 5161   3rd Qu.:0.15325   3rd Qu.:0.4098   3rd Qu.:33.25  
##  Max.   :20772   Max.   :0.33300   Max.   :0.8710   Max.   :47.00  
##   PMCiclistas       PMTPublico       PMVMotor        PeatAuto     
##  Min.   : 0.460   Min.   : 8.00   Min.   :17.00   Min.   :  5.00  
##  1st Qu.: 1.000   1st Qu.:11.75   1st Qu.:43.72   1st Qu.: 18.25  
##  Median : 2.000   Median :19.00   Median :53.00   Median : 36.50  
##  Mean   : 2.811   Mean   :21.84   Mean   :51.39   Mean   : 70.71  
##  3rd Qu.: 3.000   3rd Qu.:31.50   3rd Qu.:63.25   3rd Qu.: 67.25  
##  Max.   :14.000   Max.   :39.00   Max.   :71.00   Max.   :398.00  
##     CicAuto           V2RMSM          V2RMAuto          AutoSM      
##  Min.   :  1.00   Min.   :  1.00   Min.   :  2.00   Min.   :  0.00  
##  1st Qu.:  5.75   1st Qu.:  3.50   1st Qu.: 15.50   1st Qu.:  4.75  
##  Median : 14.00   Median :  6.00   Median : 29.50   Median :  7.50  
##  Mean   : 30.46   Mean   : 21.96   Mean   : 74.83   Mean   : 18.04  
##  3rd Qu.: 27.25   3rd Qu.: 20.00   3rd Qu.: 65.25   3rd Qu.: 21.00  
##  Max.   :340.00   Max.   :143.00   Max.   :371.00   Max.   :138.00  
##     AutoAuto      Precipitacion         Temp             PBI        
##  Min.   :  1.00   Min.   : 420.9   Min.   : 6.800   Min.   : 23867  
##  1st Qu.:  6.00   1st Qu.: 637.9   1st Qu.: 9.925   1st Qu.: 34399  
##  Median : 17.50   Median : 800.3   Median :11.100   Median : 37087  
##  Mean   : 40.17   Mean   : 753.5   Mean   :12.021   Mean   : 50517  
##  3rd Qu.: 49.25   3rd Qu.: 829.6   3rd Qu.:14.475   3rd Qu.: 48993  
##  Max.   :217.00   Max.   :1245.3   Max.   :18.200   Max.   :152178



Estudio de las covariables relacionadas con la poblacion


Vamosa empezar nuestro estudio de las covariables poniendo foco en la población de cada ciudad

boxplot(SegVial$Poblacion, main = "Poblacion", col = "blue", horizontal = T)

A priori, vemos cómo existen algunas ciudades que sobresalen en cuanto a población se refiere.

SegVial %>% ggplot(aes(x=reorder(Ciudad,Poblacion,max), y=Poblacion))+ geom_bar(stat = 'identity', col='red') +
  geom_bar(stat = 'identity', aes(x=Ciudad, y=PBI), col='blue') +
  theme(axis.text.x=element_text(angle = 90)) +
  labs(title="Población y PBI en cada ciudad")+
  xlab('Ciudades')

En este gráfico podemos ver la Población de cada ciudad en rojo y el PBI de las mismas en azul. Observamos que las 4 ciudades de las que hablábamos son: Inner London, Madrid, Rome y Paris, teniendo Inner London y París 2 de los PBI más elevados. Quizá el caso que más llama la atención es Montpellier, que presenta baja población pero muy alto PBI.

aux_plot = copy(SegVial)
aux_plot$PBI_per_capita = aux_plot$PBI / aux_plot$Poblacion
aux_plot %>% ggplot(aes(x=reorder(Ciudad,PBI_per_capita,max), y=PBI_per_capita))+ geom_bar(stat = 'identity', col='green') +
  theme(axis.text.x=element_text(angle = 90)) +
  labs(title="PBI per capita en cada ciudad")+
  xlab('Ciudades')

Vemos como al calcular el PBI per cápita Montpellier destaca y Madrid pasa a ser la ciudad con el indicador más bajo.


Analicemos ahora la Densidad Poblacional de cada ciudad

boxplot(SegVial$DenPob, main = "Densidad Poblacional", col = "blue", horizontal = T)

SegVial %>% ggplot(aes(x=reorder(Ciudad,DenPob,max), y=DenPob))+ geom_bar(stat = 'identity', col='red') +
  theme(axis.text.x=element_text(angle = 90)) +
  labs(title="Densidad Poblacional de cada ciudad")+
  xlab('Ciudades')

Observamos 4 ciudades que presentan valores de densidad poblacional mas elevadas que el resto: Paris, Barcelona, Inner London y Lyon. Sabemos que ciudades como Inner London o París tienen una elevada población, pero ciudades como Barcelona no, pero por otro lado también sabemos que la densidad poblacional es un ratio entre población sobre km2, veamos cómo son los km2 para cada ciudad.

aux_plot$km2 = aux_plot$Poblacion / aux_plot$DenPob
aux_plot %>% ggplot(aes(x=reorder(Ciudad,km2,max), y=km2))+ geom_bar(stat = 'identity', col='black') +
  theme(axis.text.x=element_text(angle = 90)) +
  labs(title="km2 de cada ciudad")+
  xlab('Ciudades')

Como decíamos antes, las ciudades con menor área van a ser más propensas a tener una concentración poblacional mayor, en este caso vemos cómo Barcelona, Lyon y Lillie tienen alta densidad poblacional no por su cuantiosa población, sino por su reducida área de extensión.



Estudio de las variables target


Este trabajo tiene como objetivo intentar darle explicación a las variables de los accidentes ocurridos, con el fin de proponer políticas públicas para disminuirlos a futuro. En este apartado vamos a investigar un poco estas features target.

Primero que nada vamos a analizar cómo son las incidencias de cada tipo de accidente en todo el dataset y luego por ciudad

color <- brewer.pal(6, "Pastel1")

pie(c(sum(SegVial$PeatAuto), sum(SegVial$CicAuto), sum(SegVial$V2RMSM),
  sum(SegVial$V2RMAuto), sum(SegVial$AutoSM), sum(SegVial$AutoAuto)), 
labels = c(" PEATON->AUTO (27.6%)"," CICLISTA->AUTO (11.9%)", " DOS RUEDAS(M)->SI MISMOS (8.6%)", 
           " DOS RUEDAS(M)->AUTO (29.2%)", " AUTO->SI MISMOS (7%)", " AUTO->AUTO (15.7%)"),
col = color, main = "DISTRIBUCION DE LOS ACCIDENTES EN EL DATASET")

Tenemos un dataset en donde la mayoría de los accidentes son de: * vehículos de 2 ruedas vs autos * autos a peatones

incidencias = copy(SegVial)
incidencias = incidencias %>% dplyr::select(Ciudad, V2RMSM, V2RMAuto, AutoAuto, AutoSM, PeatAuto, CicAuto, accidentes_viales)
incidencias$V2RMSM = round(incidencias$V2RMSM / incidencias$accidentes_viales,4)
incidencias$V2RMAuto = round(incidencias$V2RMAuto / incidencias$accidentes_viales,4)
incidencias$AutoAuto = round(incidencias$AutoAuto / incidencias$accidentes_viales,4)
incidencias$AutoSM = round(incidencias$AutoSM / incidencias$accidentes_viales,4)
incidencias$PeatAuto = round(incidencias$PeatAuto / incidencias$accidentes_viales,4)
incidencias$CicAuto = round(incidencias$CicAuto / incidencias$accidentes_viales,4)
incidencias = melt(as.data.table(incidencias))
## Warning in melt.data.table(as.data.table(incidencias)): id.vars and measure.vars
## are internally guessed when both are 'NULL'. All non-numeric/integer/
## logical type columns are considered id.vars, which in this case are columns
## [Ciudad, ...]. Consider providing at least one of 'id' or 'measure' vars in
## future.
incidencias[incidencias$variable!='accidentes_viales',] %>% ggplot(aes(x=Ciudad ,y=value, fill  = variable))+ geom_bar(stat = "identity") +
  theme(axis.text.x=element_text(angle = 90)) +
  labs(title="Incidencia del tipo de accidentes por ciudad")

De este grafico podemos sacar información interesante:


Ahora vamos a ver más en detalle los valores que toman las variables respuesta en el dataset

par(mfrow = c(3, 2))


boxplot(SegVial$PeatAuto, col = "blue",
     main = "PEATON->AUTO", horizontal = T)

boxplot(SegVial$CicAuto, col = "blue",
     main = "CICLISTA->AUTO", horizontal = T)

boxplot(SegVial$V2RMSM, col = "blue",
     main = "MOTO->SI MISMO", horizontal = T)

boxplot(SegVial$V2RMAuto, col = "blue",
     main = "MOTO->AUTO", horizontal = T)

boxplot(SegVial$AutoSM, col = "blue",
     main = "AUTO->SI MISMO", horizontal = T)

boxplot(SegVial$AutoAuto, col = "blue",
     main = "AUTO->AUTO", horizontal = T)

aux_para_plot = reshape2::melt(SegVial[,c('Ciudad','V2RMSM','V2RMAuto','AutoAuto','AutoSM','PeatAuto','CicAuto','accidentes_viales')], value.name = "accidentes")
## Using Ciudad as id variables
aux_para_plot %>% ggplot(aes(x=Ciudad ,y=accidentes, alpha=.2, color = variable))+ geom_point() +
  theme(axis.text.x=element_text(angle = 90)) +
  labs(title="Cantidad de accidentes por tipo en cada ciudad")

De estos gráfico podemos ver que: * Hay ciudades con cantidades atípicamente altas de accidentes de los diversos tipos * Inner London resulta ser la ciudad que tiene mayor cantidad de accidentes totales, accidentes PeatAuto y accidentes V2RMAuto * Roma y Madrid ocupan el 2do y 3er puesto respectivamente en cuando a accidentes totales

Recordemos del apartado de población, que estas son las 3 ciudades que mayor población tienen.

Veamos ahora con un poco más de detalle por tipo de accidente

Como dijimos, observamos que tanto la ciudad de Inner London, como Roma y Madrid aparecen en la mayoria de los casos con valores muy altos de accidentes de distintos tipos.



Estudio de variables regresoras


Areas ciclista y baja velocidad

Vamos a continuar analizando las covariables relacionadas al area destinada para cada uso

par(mfrow = c(2,1))

boxplot(SegVial$ArCiclista, main = "Area Ciclista vs Motor", col = "blue",horizontal = T)
boxplot(SegVial$ArBajaVel, main = "Area baja velocidad", col = "blue",horizontal = T)

Sobre la relacion entre el area destinada de baja velocidad, no observamos valores alejados del grupo para alguna ciudad en particular.

round(mean(SegVial$ArBajaVel), 2)
## [1] 0.3

Observamos que en promedio 7 de cada 10 km2 no tienen restricciones de velocidad en las Ciudades.

SegVial %>% 
  dplyr::select(Pais, Ciudad, ArCiclista) %>% 
  arrange(-ArCiclista) %>% head(2)
## # A tibble: 2 × 3
##   Pais   Ciudad     ArCiclista
##   <chr>  <chr>           <dbl>
## 1 France Strasbourg      0.333
## 2 France Nantes          0.327

Respecto a la relacion entre el area destinada para Ciclistas vs Vehiculos Motorizados, observamos 2 Ciudades de Francia, que presentan una cantidad de superficie mayor para el uso de la bicicleta, respecto al grupo de ciudades estudiado en general.


Parcitipación modal
aux_para_plot = reshape2::melt(SegVial[,c('Ciudad','PMPeatones','PMCiclistas','PMTPublico','PMVMotor')], value.name = "PM")
## Using Ciudad as id variables
aux_para_plot %>% ggplot(aes(x=Ciudad ,y=PM, alpha=.2, fill  = variable))+ geom_bar(stat = "identity") +
  theme(axis.text.x=element_text(angle = 90)) +
  labs(title="PM por tipo de movilidad en cada ciudad")

A simple vista podemos ver que: * Existe una clara preponderancia del transporte via vehiculo motorizado * París es la ciudad con una menor PM de vehículos a motor y una mayor incidencia de personas que se desplazan caminando. Recordemos que París es una de las ciudades con mayor cantidad de accidentes de V2RM * Barcelona también tiene poco PM de vehículos con motor, preponderando el transporte público * Leeds, Bordeaux, Roma y Sheffield tienen muy bajo PM de peatones * Bristol y Strasbourg tienen niveles elevados de ciclistas, respecto a las demas ciudades


Variables meteorológicas
par(mfrow = c(2,1))

boxplot(SegVial$Temp, main = "TEMPERATURA MEDIA [C]", col = "blue", horizontal = T)
boxplot(SegVial$Precipitacion, main = "PRECIPITACION MEDIA [mm]", col = "blue", horizontal = T)

SegVial %>% 
  dplyr::select(Pais, Ciudad, Precipitacion) %>% 
  arrange(-Precipitacion) %>% head(1)
## # A tibble: 1 × 3
##   Pais           Ciudad  Precipitacion
##   <chr>          <chr>           <dbl>
## 1 United Kingdom Glasgow         1245.

Observamos un valor elevado de mm de agua para la ciudad de Glasgow, ¿Esto impacta en la cantidad de ciclistas?

par(mfrow = c(2, 2))


plot(y = SegVial$PMCiclistas, x = SegVial$Precipitacion, col = "blue",
     main = "modal CICLISTAS vs mm PRECIPITADOS", 
     xlab = "mm PRECIPITADOS", ylab = "CICLISTAS [%]")

plot(y = SegVial$PMPeatones, x = SegVial$Precipitacion, col = "blue",
     main = "modal PEATONES vs mm PRECIPITADOS", 
     xlab = "mm PRECIPITADOS", ylab = "PEATONES [%]")


plot(y = SegVial$PMCiclistas, x = SegVial$Temp, col = "blue",
     main = "modal CICLISTAS vs Temp Media")

plot(y = SegVial$PMPeatones, x = SegVial$Temp, col = "blue",
     main = "modal PEATONES vs Temp Media")

cor(SegVial$PMCiclistas, SegVial[, c(18:19)], method = "spearman")
##      Precipitacion        Temp
## [1,]   -0.04857822 -0.08844308

No se observa una correlacion entre una ciudad con muchas lluvias o con mucha temperatura y el optar o no por el transporte en bicicleta.

Revisemos lo mismo con peatones

cor(SegVial$PMPeatones, SegVial[, c(18:19)], method = "spearman")
##      Precipitacion     Temp
## [1,]    -0.5121323 0.345434

Aqui si se observa una correlacion un poco más fuerte. Para las ciudades que tienen una media de temperatura mas elevada, la gente utiliza mas el transporte a pie, lo cual tiene sentido. A su vez, tambien observamos que ciudades que presentan mas cantidad de lluvias, ven reducido el uso del transporte a pie.

cor(SegVial$PMTPublico, SegVial[, c(18:19)], method = "spearman")
##      Precipitacion       Temp
## [1,]     0.1179801 -0.3302114
cor(SegVial$PMVMotor, SegVial[, c(18:19)], method = "spearman")
##      Precipitacion        Temp
## [1,]     0.3628454 -0.03156291

Observamos que ciudades con mayor cantidad de precipitaciones, y temperaturas medias mas bajas, tienen un aumento del transporte motorizado.



Continuamos el analisis respondiendo algunas preguntas típicas de los accidentes viales:

¿Mayor densidad de poblacion favorece la cantidad de accidentes en general?

par(mfrow = c(2,3))
plot(SegVial$DenPob, SegVial$AutoAuto, main='DenPob vs accidentes Auto-Auto')
plot(SegVial$DenPob, SegVial$AutoSM, main='DenPob vs accidentes AutoSM')
plot(SegVial$DenPob, SegVial$V2RMAuto, main='DenPob vs accidentes V2RM-Auto')
plot(SegVial$DenPob, SegVial$V2RMSM, main='DenPob vs accidentes V2RMSM')
plot(SegVial$DenPob, SegVial$PeatAuto, main='DenPob vs accidentes Peat-Auto')
plot(SegVial$DenPob, SegVial$CicAuto, main='DenPob vs accidentes Cic-Auto')

Por sentido común, esperaríamos que poblaciones con mayor densidad sean propensas a más accidentes, de cualquier tipo, ya que la concentración de gente en zonas más pequeñas deberia impactar en el ordenamiento del tránsito. Sin embargo, los plots no indican esto, por lo que DenPob no pareciera ser una variable muy indicadora del nivel de accidentes.

¿Mayor cantidad de poblacion favorece a la cantidad de accidentes?

par(mfrow = c(2,3))
plot(SegVial$Poblacion, SegVial$AutoAuto, main='Poblacion vs accidentes Auto-Auto')
plot(SegVial$Poblacion, SegVial$AutoSM, main='Poblacion vs accidentes AutoSM')
plot(SegVial$Poblacion, SegVial$V2RMAuto, main='Poblacion vs accidentes V2RM-Auto')
plot(SegVial$Poblacion, SegVial$V2RMSM, main='Poblacion vs accidentes V2RMSM')
plot(SegVial$Poblacion, SegVial$PeatAuto, main='Poblacion vs accidentes Peat-Auto')
plot(SegVial$Poblacion, SegVial$CicAuto, main='Poblacion vs accidentes Cic-Auto')

t(cor(SegVial$Poblacion, SegVial[, c(12:17)], 
                         method = "spearman"))
##               [,1]
## PeatAuto 0.7598957
## CicAuto  0.5587469
## V2RMSM   0.6363902
## V2RMAuto 0.7401609
## AutoSM   0.6470724
## AutoAuto 0.6930785

Observamos que la población tiene relación aproximadamente lineal respecto a los distintos tipos de accidentes. No sabemos si las variables de accidente son realmente índices o totales de accidentes. Si fueran lo segundo, tiene sentido el efecto - donde hay más gente, podés reportar más casos-. Sin embargo, cabe destacar como vimos en los gráficos previos, que esa gente esté más o menos distante, impacta poco.

nota: Cabe aclarar que se evidencia visualmente la presencia de valores atipicamente altos, que fueron analizados previamente en la seccion “ANALISIS DE LOS TARGETS”.

Vamos a analizar si alguna covariable presenta una correlacion lineal contra la variable target_combinacion (sumatoria de todos los accidentes por ciudad). Si fuera asi, esto indicaria que tiene algo que explicar con respecto a alguno de los targets.

df_pairs_1 = SegVial[,-c(1:3, 5, 12:20)]

ggpairs(df_pairs_1, progress = F)

df_pairs_2 = SegVial[,c(18:21)]

ggpairs(df_pairs_2, progress = F)

Observamos una gran correlacion entre nuestro target_combination y poblacion, asi como tambien con PMTPublico y a continuacion lo siguen ArCiclista, Precipitacion y PMCiclistas. Estas covariables probablemente puedan explicarnos bastante de nuestros targets al momento de armar los modelos.



SegVial$PeatAuto_std = (SegVial$PeatAuto / SegVial$Poblacion) * (10 ^ 6)
SegVial$CicAuto_std  = (SegVial$CicAuto  / SegVial$Poblacion) * (10 ^ 6)
SegVial$V2RMSM_std   = (SegVial$V2RMSM   / SegVial$Poblacion) * (10 ^ 6)
SegVial$V2RMAuto_std = (SegVial$V2RMAuto / SegVial$Poblacion) * (10 ^ 6)
SegVial$AutoSM_std   = (SegVial$AutoSM   / SegVial$Poblacion) * (10 ^ 6)
SegVial$AutoAuto_std = (SegVial$AutoAuto / SegVial$Poblacion) * (10 ^ 6)
SegVial$accidentes_viales_std = (SegVial$accidentes_viales / SegVial$Poblacion) * (10 ^ 6)


SegVial[,-c(1:3, 12:17, 21:26)] = as.data.frame(scale(SegVial[,-c(1:3, 12:17, 21:26)]))

head(SegVial)
## # A tibble: 6 × 28
##   `Codigo Pais` Pais   Ciudad  Poblacion  DenPob ArCiclista ArBajaVel PMPeatones
##   <chr>         <chr>  <chr>       <dbl>   <dbl>      <dbl>     <dbl>      <dbl>
## 1 ES            Spain  Barcel…     0.654  2.23       0.404     1.90        0.662
## 2 ES            Spain  Madrid      2.28  -0.0320    -0.605    -0.924       0.827
## 3 FR            France Bordea…    -0.734 -0.0802    -0.131    -0.232      -0.246
## 4 FR            France Lille      -0.756  0.255     -0.228     1.34        0.662
## 5 FR            France Lyon       -0.468  1.12      -0.0942   -0.0533      0.827
## 6 FR            France Marsei…    -0.116 -0.406     -0.884    -0.713       0.827
## # … with 20 more variables: PMCiclistas <dbl>, PMTPublico <dbl>,
## #   PMVMotor <dbl>, PeatAuto <dbl>, CicAuto <dbl>, V2RMSM <dbl>,
## #   V2RMAuto <dbl>, AutoSM <dbl>, AutoAuto <dbl>, Precipitacion <dbl>,
## #   Temp <dbl>, PBI <dbl>, accidentes_viales <dbl>, PeatAuto_std <dbl>,
## #   CicAuto_std <dbl>, V2RMSM_std <dbl>, V2RMAuto_std <dbl>, AutoSM_std <dbl>,
## #   AutoAuto_std <dbl>, accidentes_viales_std <dbl>

Como ya dijimos, el objetivo del presente trabajo es encontrar caracteristicas urbanas para encontrar las vulnerabilidades de los usuarios viales, y de esta manera proponer politicas publicas eficaces para disminuir la fatalidad y cantidad de accidentes. Dicho esto, dado las características del dataset, se evidencia que las covariables Pais, Codigo Pais y Ciudad, no nos darian informacion util para este analisis, y por tanto optamos por descartarlas.

df_peatauto = SegVial[, -c(1:3, 12:17, 21, 23:28)]
df_cicauto  = SegVial[, -c(1:3, 12:17, 21, 22, 24:28)]
df_v2rmsm   = SegVial[, -c(1:3, 12:17, 21, 23, 25:28)]
df_v2rmauto = SegVial[, -c(1:3 ,12:17, 21:24, 26:28)]
df_autosm   = SegVial[, -c(1:3 ,12:17, 21:25, 27:28)]
df_autoauto = SegVial[, -c(1:3 ,12:17, 21:26, 28)]
df_accidentesviales = SegVial[, -c(1:3 ,12:17, 21:27)]




TARGET PEATAUTO


modelo_peatauto_ls = lm(PeatAuto_std ~., df_peatauto)

summary(modelo_peatauto_ls)
## 
## Call:
## lm(formula = PeatAuto_std ~ ., data = df_peatauto)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -45.947  -9.491  -0.023  13.550  31.458 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    67.2207     5.3963  12.457 3.19e-08 ***
## Poblacion       9.6839     8.8687   1.092   0.2963    
## DenPob          1.2058     9.7536   0.124   0.9037    
## ArCiclista    -11.9114     7.2237  -1.649   0.1251    
## ArBajaVel     -11.3210     8.2395  -1.374   0.1946    
## PMPeatones    436.1806   220.4523   1.979   0.0713 .  
## PMCiclistas   102.0832    52.6049   1.941   0.0762 .  
## PMTPublico    377.1615   184.6464   2.043   0.0637 .  
## PMVMotor      514.9586   250.2670   2.058   0.0620 .  
## Precipitacion  -0.2775     8.1269  -0.034   0.9733    
## Temp           -8.7025     8.2750  -1.052   0.3137    
## PBI            -1.3527     6.7087  -0.202   0.8436    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 26.44 on 12 degrees of freedom
## Multiple R-squared:  0.7489, Adjusted R-squared:  0.5187 
## F-statistic: 3.254 on 11 and 12 DF,  p-value: 0.02693

Observamos un R2adj del orden del 50%, y un pvalue global, no tan bueno como nos gustaria, pero significativo. Tambien notamos que varios estimadores no resultan significativamente distintos de 0, esto puede deberse a problemas de alta colinealidad entre covariables. Vamos a realizar un analisis de la colinealidad entre las covariables.


Vamos a realizar un analisis VIF


nota: Aplicamos una transformacion logaritmica para una visualizacion mas amigable, los valores por sobre la linea roja, indican correlacion elevada.

barplot(log(vif(modelo_peatauto_ls))[vif(modelo_peatauto_ls) > 5],
          main = "log VIF DE LAS COVARIABLES", col = "blue")

abline(h = log(5), lwd = 2, col = "red")

(Cabe destacar en el gráfico se muestra el logaritmo de los valores del VIF, para facilitar la visualización)

Se observan problemas de colinealidad entre las participaciones modales de las formas de transporte: PMPeatones, PMCiclistas, PMTPublic y PMVmotor. Esto tiene logica desde la construccion de las variables modales.


Vamos a trabajar para intentar solucionar los problemas de multicolinealidad hallados.


Un primer approach podemos hacerlo a partir de metodos de seleccion de covariables como el criterio stepwise.

METODO STEPWISE TARGET PEATAUTO

modelo_peatauto_wise=step(modelo_peatauto_ls, direction = c("both"), trace = F)

summary(modelo_peatauto_wise)
## 
## Call:
## lm(formula = PeatAuto_std ~ Poblacion + ArCiclista + ArBajaVel + 
##     PMPeatones + PMCiclistas + PMTPublico + PMVMotor + Temp, 
##     data = df_peatauto)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -47.897  -9.821   0.269  14.257  33.062 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   67.221      4.836  13.899 5.66e-10 ***
## Poblacion      9.330      6.850   1.362   0.1933    
## ArCiclista   -12.101      6.240  -1.939   0.0715 .  
## ArBajaVel    -10.705      6.491  -1.649   0.1199    
## PMPeatones   424.800    183.052   2.321   0.0348 *  
## PMCiclistas   99.486     43.991   2.262   0.0390 *  
## PMTPublico   368.083    153.363   2.400   0.0298 *  
## PMVMotor     501.859    208.362   2.409   0.0293 *  
## Temp          -7.922      6.092  -1.300   0.2131    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 23.69 on 15 degrees of freedom
## Multiple R-squared:  0.7479, Adjusted R-squared:  0.6135 
## F-statistic: 5.563 on 8 and 15 DF,  p-value: 0.002166

Pasamos de tener 11 a 8 covariables en nuestro modelo, pero observamos que el metodo opto por quedarse con las covariables de participacion modal. Tambien, observamos una mejoria en el R2adj, asi como tambien en el RSE. Tambien logramos aumentar la cantidad de estimadores estadisticamente distintos de cero.

barplot(log(vif(modelo_peatauto_wise))[vif(modelo_peatauto_wise) > 5],
          main = "log VIF DE LAS COVARIABLES", col = "blue")

abline(h = log(5), lwd = 2, col = "red")

Seguimos manteniendo el problema de multicolinealidad, vamos a probar otro enfoque.

Vamos a aplicar un metodo de regularizacion, para quitar covariables, a ver si logramos mejorar nuestro modelo.


METODO DE REGULARIZACION LASSO TARGET PEATAUTO

xx = model.matrix(PeatAuto_std ~., data = df_peatauto)[, -1]

modelo_peatauto_lassocv <- cv.glmnet(xx, y = df_peatauto$PeatAuto_std, 
                                         alpha = 1,
                                         nfolds = 5)
modelo_peatauto_lassocv_min <- glmnet(xx, y = df_peatauto$PeatAuto_std, 
                                         alpha = 1,
                                         lambda = modelo_peatauto_lassocv$lambda.min)

modelo_peatauto_lassocv_1sd <- glmnet(xx, y = df_peatauto$PeatAuto_std, 
                                         alpha = 1,
                                         lambda = modelo_peatauto_lassocv$lambda.1se)

cbind(coef(modelo_peatauto_lassocv_min),
      coef(modelo_peatauto_lassocv_1sd))
## 12 x 2 sparse Matrix of class "dgCMatrix"
##                      s0        s0
## (Intercept)    67.22071 67.220709
## Poblacion       .        .       
## DenPob          .        .       
## ArCiclista    -11.73146 -5.340420
## ArBajaVel       .        .       
## PMPeatones    -13.59697 -7.205842
## PMCiclistas     .        .       
## PMTPublico      .        .       
## PMVMotor        .        .       
## Precipitacion   .        .       
## Temp            .        .       
## PBI             .        .

Ambos criterios de penalidad deciden conservar las mismas covariables para la regresion

modelo_peatauto_ls_2=lm(PeatAuto_std ~ ArCiclista + PMPeatones, data = df_peatauto)

summary(modelo_peatauto_ls_2)
## 
## Call:
## lm(formula = PeatAuto_std ~ ArCiclista + PMPeatones, data = df_peatauto)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -37.304 -22.595   2.063  14.025  43.113 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   67.221      5.192  12.948 1.77e-11 ***
## ArCiclista   -16.609      5.781  -2.873  0.00911 ** 
## PMPeatones   -18.475      5.781  -3.196  0.00435 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 25.43 on 21 degrees of freedom
## Multiple R-squared:  0.5933, Adjusted R-squared:  0.5545 
## F-statistic: 15.32 on 2 and 21 DF,  p-value: 7.899e-05

Observamos pvalue global e individuales mucho mas significativos. Nuestro R2adj se ubica en un 55% para el target PeatAuto_std.

par(mfrow=c(2,2))
plot(modelo_peatauto_ls_2)

Revisamos si los residuos poseen algun tipo de estructura o si se encuentra algun outlier, pero todo parece estar en orden.

control = lmrobdet.control(bb = 0.5, efficiency = 0.85, family = "bisquare")

modelo_MM_peatauto= lmrobdetMM(PeatAuto_std ~ ArCiclista + 
                              PMPeatones, data = df_peatauto, control = control)

summary(modelo_MM_peatauto)
## 
## Call:
## lmrobdetMM(formula = PeatAuto_std ~ ArCiclista + PMPeatones, data = df_peatauto, 
##     control = control)
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -37.507 -20.799   2.045  15.554  43.535 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   66.761      6.060  11.018 3.46e-10 ***
## ArCiclista   -16.000      6.561  -2.438   0.0237 *  
## PMPeatones   -17.263      6.565  -2.630   0.0157 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Robust residual standard error: 30.81 
## Multiple R-squared:  0.6318, Adjusted R-squared:  0.5967 
## Convergence in 11 IRWLS iterations
plot(modelo_MM_peatauto$residuals, modelo_MM_peatauto$rweights,
     main = "ESTUDIO DE OUTLIERS", xlab = "RESIDUOS", ylab = "RWEIGHTS", col = "blue")

No se observan outliers, pero conseguimos un mejor R2adj que con el metodo de LS.



TARGET CICAUTO


modelo_cicauto_ls = lm(CicAuto_std ~., data = df_cicauto)

summary(modelo_cicauto_ls)
## 
## Call:
## lm(formula = CicAuto_std ~ ., data = df_cicauto)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -32.807  -9.711   4.135   9.997  24.201 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    29.2005     4.4369   6.581 2.61e-05 ***
## Poblacion       3.7392     7.2920   0.513    0.617    
## DenPob          0.1312     8.0196   0.016    0.987    
## ArCiclista     -3.2262     5.9395  -0.543    0.597    
## ArBajaVel       2.1050     6.7747   0.311    0.761    
## PMPeatones    233.0700   181.2602   1.286    0.223    
## PMCiclistas    63.8725    43.2528   1.477    0.166    
## PMTPublico    208.2031   151.8199   1.371    0.195    
## PMVMotor      278.5934   205.7745   1.354    0.201    
## Precipitacion  -2.0798     6.6821  -0.311    0.761    
## Temp           -9.3761     6.8038  -1.378    0.193    
## PBI            -0.0853     5.5161  -0.015    0.988    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 21.74 on 12 degrees of freedom
## Multiple R-squared:  0.6285, Adjusted R-squared:  0.288 
## F-statistic: 1.846 on 11 and 12 DF,  p-value: 0.1535
barplot(log(vif(modelo_cicauto_ls))[vif(modelo_cicauto_ls) > 5],
          main = "log VIF DE LAS COVARIABLES", col = "blue")

abline(h = log(5), lwd = 2, col = "red")


METODO STEPWISE TARGET CICAUTO


modelo_cicauto_wise=step(modelo_cicauto_ls, direction = c("both"), trace = F)

summary(modelo_cicauto_wise)
## 
## Call:
## lm(formula = CicAuto_std ~ PMPeatones + PMCiclistas + PMTPublico + 
##     PMVMotor + Temp, data = df_cicauto)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -33.058 -11.882   2.272   7.486  36.046 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   29.201      3.832   7.621 4.86e-07 ***
## PMPeatones   244.590    125.882   1.943   0.0678 .  
## PMCiclistas   65.834     30.005   2.194   0.0416 *  
## PMTPublico   220.626    106.099   2.079   0.0521 .  
## PMVMotor     291.050    144.061   2.020   0.0585 .  
## Temp          -6.907      4.424  -1.561   0.1359    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 18.77 on 18 degrees of freedom
## Multiple R-squared:  0.5845, Adjusted R-squared:  0.469 
## F-statistic: 5.064 on 5 and 18 DF,  p-value: 0.004519


METODO DE REGULARIZACION LASSO TARGET CICAUTO


xx = model.matrix(CicAuto_std ~., data = df_cicauto)[, -1]

modelo_cicauto_lassocv <- cv.glmnet(xx, y = df_cicauto$CicAuto_std, 
                                         alpha = 1,
                                         nfolds = 5)
modelo_cicauto_lassocv_min <- glmnet(xx, y = df_cicauto$CicAuto_std, 
                                         alpha = 1,
                                         lambda = modelo_cicauto_lassocv$lambda.min)

modelo_cicauto_lassocv_1sd <- glmnet(xx, y = df_cicauto$CicAuto_std, 
                                         alpha = 1,
                                         lambda = modelo_cicauto_lassocv$lambda.1se)

cbind(coef(modelo_cicauto_lassocv_min),
      coef(modelo_cicauto_lassocv_1sd))
## 12 x 2 sparse Matrix of class "dgCMatrix"
##                       s0       s0
## (Intercept)   29.2005120 29.20051
## Poblacion      .          0.00000
## DenPob         .          .      
## ArCiclista     .          .      
## ArBajaVel      .          .      
## PMPeatones    -4.6539351  .      
## PMCiclistas    .          .      
## PMTPublico     0.1524345  .      
## PMVMotor       .          .      
## Precipitacion  .          .      
## Temp          -5.7587640  .      
## PBI            .          .

Observamos que el criterio 1sd de lambda, directamente nos deja sin covariables y que el criterio min, nos sugiere que dejemos tanto PMPeatones como PMTPublico, dado que ambas tienen colinealidad como marcamos previamente, optamos por dejar solo PMPeatones

modelo_cicauto_ls_2=lm(CicAuto_std ~ PMPeatones + Temp, data = df_cicauto)

summary(modelo_cicauto_ls_2)
## 
## Call:
## lm(formula = CicAuto_std ~ PMPeatones + Temp, data = df_cicauto)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -34.042  -4.402  -0.553   5.602  63.282 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   29.201      4.210   6.935  7.5e-07 ***
## PMPeatones    -9.680      4.529  -2.137   0.0445 *  
## Temp         -10.776      4.529  -2.379   0.0269 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 20.63 on 21 degrees of freedom
## Multiple R-squared:  0.4146, Adjusted R-squared:  0.3589 
## F-statistic: 7.437 on 2 and 21 DF,  p-value: 0.003616

No nos queda un modelo muy bueno, por ahora nos quedariamos con modelo_cicauto_wise

par(mfrow=c(2,2))
plot(modelo_cicauto_wise)

Se observan registros a tener en cuenta con alto leverage y residuo cercano a 3

modelo_MM_cicauto = lmrobdetMM(CicAuto_std ~ PMPeatones + Temp, data = df_cicauto, control = control)

summary(modelo_MM_cicauto)
## 
## Call:
## lmrobdetMM(formula = CicAuto_std ~ PMPeatones + Temp, data = df_cicauto, 
##     control = control)
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -35.7783  -7.2515  -0.4663   4.4901  61.3931 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   31.090      2.961  10.500  8.2e-10 ***
## PMPeatones   -12.388      3.498  -3.542  0.00193 ** 
## Temp         -11.656      3.556  -3.278  0.00359 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Robust residual standard error: 12.02 
## Multiple R-squared:  0.4775, Adjusted R-squared:  0.4277 
## Convergence in 68 IRWLS iterations
plot(modelo_MM_cicauto$residuals, modelo_MM_cicauto$rweights,
     main = "ESTUDIO DE OUTLIERS", xlab = "RESIDUOS", ylab = "RWEIGHTS", col = "blue")

Encontramos 1 outlier, al que la regresion robusta le asigna un rweight de 0, y que tiene un residuo muy elevado.

sort(modelo_MM_cicauto$rweights)[1]
## 20 
##  0

Dropeamos esta ciudad y creamos un modelo OLS:

modelo_cicauto_ls_3 = lm(CicAuto_std ~ PMPeatones + Temp, data = df_cicauto[-c(20),])

summary(modelo_cicauto_ls_3)
## 
## Call:
## lm(formula = CicAuto_std ~ PMPeatones + Temp, data = df_cicauto[-c(20), 
##     ])
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -29.246  -2.486   1.355   6.310  26.473 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   26.437      3.210   8.237  7.4e-08 ***
## PMPeatones   -10.205      3.381  -3.018  0.00680 ** 
## Temp          -9.870      3.386  -2.915  0.00856 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 15.39 on 20 degrees of freedom
## Multiple R-squared:  0.5622, Adjusted R-squared:  0.5184 
## F-statistic: 12.84 on 2 and 20 DF,  p-value: 0.0002587



TARGET V2RMSM


modelo_v2rmsm_ls = lm(V2RMSM_std ~., data = df_v2rmsm)

summary(modelo_v2rmsm_ls)
## 
## Call:
## lm(formula = V2RMSM_std ~ ., data = df_v2rmsm)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -11.3564  -5.9630   0.5517   4.1226  16.6726 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)  
## (Intercept)    22.78329    7.85880   2.899   0.0145 *
## Poblacion       5.79059    3.62827   1.596   0.1388  
## DenPob         -0.61035    3.80811  -0.160   0.8756  
## ArCiclista     -0.71819    3.12160  -0.230   0.8223  
## ArBajaVel      -7.88648    3.45857  -2.280   0.0435 *
## PMPeatones    145.09200   99.05874   1.465   0.1710  
## PMCiclistas    33.30006   23.52679   1.415   0.1846  
## PMTPublico    122.81272   83.63822   1.468   0.1700  
## PMVMotor      167.35331  113.57767   1.473   0.1687  
## Precipitacion  -3.69178    3.17114  -1.164   0.2690  
## Temp            4.99689    3.37427   1.481   0.1667  
## PBI             0.18098    2.62207   0.069   0.9462  
## PeatAuto_std   -0.09294    0.11264  -0.825   0.4268  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 10.32 on 11 degrees of freedom
## Multiple R-squared:  0.7151, Adjusted R-squared:  0.4043 
## F-statistic: 2.301 on 12 and 11 DF,  p-value: 0.08913
barplot(log(vif(modelo_v2rmsm_ls))[vif(modelo_v2rmsm_ls) > 5],
          main = "log VIF DE LAS COVARIABLES", col = "blue")

abline(h = log(5), lwd = 2, col = "red")


METODO STEPWISE TARGET V2RMSM


modelo_v2rmsm_wise=step(modelo_v2rmsm_ls, direction = c("both"), trace = F)

summary(modelo_v2rmsm_wise)
## 
## Call:
## lm(formula = V2RMSM_std ~ Poblacion + ArBajaVel + PMPeatones + 
##     PMCiclistas + PMTPublico + PMVMotor + Precipitacion + Temp, 
##     data = df_v2rmsm)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -11.6412  -6.6659  -0.4107   4.1672  16.8782 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     16.536      1.863   8.876 2.34e-07 ***
## Poblacion        4.818      2.890   1.667   0.1162    
## ArBajaVel       -7.089      2.554  -2.775   0.0142 *  
## PMPeatones     105.585     71.474   1.477   0.1603    
## PMCiclistas     24.200     17.200   1.407   0.1798    
## PMTPublico      88.430     60.056   1.472   0.1616    
## PMVMotor       120.774     81.395   1.484   0.1586    
## Precipitacion   -3.854      2.721  -1.416   0.1771    
## Temp             5.468      2.446   2.235   0.0410 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.127 on 15 degrees of freedom
## Multiple R-squared:  0.6959, Adjusted R-squared:  0.5337 
## F-statistic:  4.29 on 8 and 15 DF,  p-value: 0.007388


METODO DE REGULARIZACION LASSO TARGET V2RMSM


xx = model.matrix(V2RMSM_std ~., data = df_v2rmsm)[, -1]

modelo_v2rmsm_lassocv <- cv.glmnet(xx, y = df_v2rmsm$V2RMSM_std, 
                                         alpha = 1,
                                         nfolds = 5)
modelo_v2rmsm_lassocv_min <- glmnet(xx, y = df_v2rmsm$V2RMSM_std, 
                                         alpha = 1,
                                         lambda = modelo_v2rmsm_lassocv$lambda.min)

modelo_v2rmsm_lassocv_1sd <- glmnet(xx, y = df_v2rmsm$V2RMSM_std, 
                                         alpha = 1,
                                         lambda = modelo_v2rmsm_lassocv$lambda.1se)

cbind(coef(modelo_v2rmsm_lassocv_min),
      coef(modelo_v2rmsm_lassocv_1sd))
## 13 x 2 sparse Matrix of class "dgCMatrix"
##                       s0         s0
## (Intercept)   16.5359592 16.5359592
## Poblacion      3.1086528  0.0722265
## DenPob         .          .        
## ArCiclista     .          .        
## ArBajaVel     -3.5556674  .        
## PMPeatones     .          .        
## PMCiclistas   -0.3550315  .        
## PMTPublico     .          .        
## PMVMotor       .          .        
## Precipitacion -1.9226689  .        
## Temp           5.0761962  3.0827599
## PBI            .          .        
## PeatAuto_std   .          .

Vamos a adoptar el criterio min en este caso

modelo_v2rmsm_ls_2=lm(V2RMSM_std ~ Poblacion + ArBajaVel + PMCiclistas+
                        Precipitacion + Temp, data = df_v2rmsm)

summary(modelo_v2rmsm_ls_2)
## 
## Call:
## lm(formula = V2RMSM_std ~ Poblacion + ArBajaVel + PMCiclistas + 
##     Precipitacion + Temp, data = df_v2rmsm)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -14.5083  -6.7125   0.1061   3.5711  20.1013 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     16.536      1.825   9.061 3.98e-08 ***
## Poblacion        4.061      2.158   1.882   0.0761 .  
## ArBajaVel       -5.375      2.041  -2.634   0.0169 *  
## PMCiclistas     -1.039      2.034  -0.511   0.6156    
## Precipitacion   -3.227      2.439  -1.323   0.2023    
## Temp             5.461      2.248   2.429   0.0258 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.94 on 18 degrees of freedom
## Multiple R-squared:  0.6498, Adjusted R-squared:  0.5525 
## F-statistic:  6.68 on 5 and 18 DF,  p-value: 0.001108
par(mfrow=c(2,2))
plot(modelo_v2rmsm_ls_2)

Hay algunas ciudades con alta distancia de cook y residuo superior a modulo de 2.

modelo_MM_v2rmsm = lmrobdetMM(V2RMSM_std ~ Poblacion + ArBajaVel + PMCiclistas + 
    Precipitacion + Temp, data = df_v2rmsm, control = control)

summary(modelo_MM_v2rmsm)
## 
## Call:
## lmrobdetMM(formula = V2RMSM_std ~ Poblacion + ArBajaVel + PMCiclistas + Precipitacion + 
##     Temp, data = df_v2rmsm, control = control)
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7.7589 -3.4036 -0.4829  4.6299 26.8394 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    14.5173     1.6220   8.950 4.78e-08 ***
## Poblacion       3.3130     2.0732   1.598   0.1275    
## ArBajaVel      -3.3675     1.7910  -1.880   0.0764 .  
## PMCiclistas    -0.8955     1.6934  -0.529   0.6034    
## Precipitacion  -2.5392     2.2508  -1.128   0.2741    
## Temp            3.6217     1.9652   1.843   0.0819 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Robust residual standard error: 8.773 
## Multiple R-squared:  0.5196, Adjusted R-squared:  0.3861 
## Convergence in 37 IRWLS iterations
plot(modelo_MM_v2rmsm$residuals, modelo_MM_v2rmsm$rweights,
     main = "ESTUDIO DE OUTLIERS", xlab = "RESIDUOS", ylab = "RWEIGHTS", col = "blue")

sort(modelo_MM_v2rmsm$rweights)[1]
##           6 
## 0.001668902
SegVial[c(6) ,"Ciudad"]
## # A tibble: 1 × 1
##   Ciudad   
##   <chr>    
## 1 Marseille

Encontramos que la ciudad de Marseille, tiene un residuo muy alto y peso bajo segun nuestro modelo robusto. Vamos a dropearla y crear otro modelo LS, sin esta observacion.

modelo_v2rmsm_ls_3 = lm(V2RMSM_std ~ Poblacion + ArBajaVel + PMCiclistas + 
                      Precipitacion + Temp, data = df_v2rmsm[-c(6),])


summary(modelo_v2rmsm_ls_3)
## 
## Call:
## lm(formula = V2RMSM_std ~ Poblacion + ArBajaVel + PMCiclistas + 
##     Precipitacion + Temp, data = df_v2rmsm[-c(6), ])
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -9.5434 -5.4103  0.0495  3.8454 13.9014 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    15.4401     1.5337  10.067  1.4e-08 ***
## Poblacion       5.2117     1.8033   2.890   0.0102 *  
## ArBajaVel      -4.2813     1.7061  -2.509   0.0225 *  
## PMCiclistas    -0.4932     1.6738  -0.295   0.7718    
## Precipitacion  -1.0105     2.1171  -0.477   0.6392    
## Temp            4.9617     1.8466   2.687   0.0156 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.317 on 17 degrees of freedom
## Multiple R-squared:  0.6734, Adjusted R-squared:  0.5774 
## F-statistic: 7.011 on 5 and 17 DF,  p-value: 0.001008

Vamos a ver si podemos encontrar alguna transformacion de las covariables que mejore aun mas nuestro modelo.

# 
df_v2rmsm$PMCiclistas = df_v2rmsm$PMCiclistas + 1
df_v2rmsm$Precipitacion = df_v2rmsm$Precipitacion + 2.1
df_v2rmsm$Temp = df_v2rmsm$Temp + 2
# 
# boxTidwell(V2RMSM_std ~ Precipitacion + PMCiclistas + Temp, data = df_v2rmsm[-c(6), ])
df_v2rmsm$PMCiclistas2 = 1/sqrt(df_v2rmsm$PMCiclistas)

df_v2rmsm$Temp = sqrt(df_v2rmsm$Temp)

df_v2rmsm$Precipitacion2 = df_v2rmsm$Precipitacion^2

modelo_v2rmsm_ls_3 = lm(V2RMSM_std ~ Poblacion + ArBajaVel + 
                     PMCiclistas2 + Precipitacion2 + Temp , data = df_v2rmsm[-c(6),])


summary(modelo_v2rmsm_ls_3)
## 
## Call:
## lm(formula = V2RMSM_std ~ Poblacion + ArBajaVel + PMCiclistas2 + 
##     Precipitacion2 + Temp, data = df_v2rmsm[-c(6), ])
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -9.6600 -5.3257 -0.1406  4.0773 13.1634 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)   
## (Intercept)      2.1745     7.9853   0.272  0.78866   
## Poblacion        5.7118     1.7797   3.209  0.00514 **
## ArBajaVel       -5.3317     1.8920  -2.818  0.01185 * 
## PMCiclistas2    -2.3264     4.8725  -0.477  0.63912   
## Precipitacion2  -0.3877     0.3343  -1.160  0.26225   
## Temp            13.3211     4.4019   3.026  0.00762 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.943 on 17 degrees of freedom
## Multiple R-squared:  0.7059, Adjusted R-squared:  0.6195 
## F-statistic: 8.163 on 5 and 17 DF,  p-value: 0.000439
modelo_v2rmsm_wise2 = step(modelo_v2rmsm_ls_3, direction = c("both"), trace = F)

summary(modelo_v2rmsm_wise2)
## 
## Call:
## lm(formula = V2RMSM_std ~ Poblacion + ArBajaVel + Precipitacion2 + 
##     Temp, data = df_v2rmsm[-c(6), ])
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -9.7252 -5.1752  0.3439  3.2655 13.7347 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)   
## (Intercept)      0.1084     6.5656   0.017  0.98701   
## Poblacion        5.3045     1.5281   3.471  0.00272 **
## ArBajaVel       -4.8303     1.5396  -3.137  0.00569 **
## Precipitacion2  -0.4162     0.3218  -1.293  0.21220   
## Temp            12.8700     4.2061   3.060  0.00675 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.792 on 18 degrees of freedom
## Multiple R-squared:  0.702,  Adjusted R-squared:  0.6358 
## F-statistic:  10.6 on 4 and 18 DF,  p-value: 0.0001356



TARGET V2RMAUTO


modelo_v2rmauto_ls = lm(V2RMAuto_std ~., data = df_v2rmauto)

summary(modelo_v2rmauto_ls)
## 
## Call:
## lm(formula = V2RMAuto_std ~ ., data = df_v2rmauto)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -38.083 -14.225  -0.007  12.132  28.776 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     62.263      4.931  12.627 2.74e-08 ***
## Poblacion       13.992      8.104   1.727   0.1099    
## DenPob          -5.643      8.913  -0.633   0.5385    
## ArCiclista     -13.765      6.601  -2.085   0.0591 .  
## ArBajaVel      -12.994      7.529  -1.726   0.1100    
## PMPeatones     514.017    201.448   2.552   0.0254 *  
## PMCiclistas    121.364     48.070   2.525   0.0267 *  
## PMTPublico     418.179    168.728   2.478   0.0290 *  
## PMVMotor       578.256    228.692   2.529   0.0265 *  
## Precipitacion   -7.792      7.426  -1.049   0.3148    
## Temp            20.817      7.562   2.753   0.0175 *  
## PBI             -2.411      6.130  -0.393   0.7010    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 24.16 on 12 degrees of freedom
## Multiple R-squared:  0.8031, Adjusted R-squared:  0.6227 
## F-statistic: 4.451 on 11 and 12 DF,  p-value: 0.008055
modelo_v2rmauto_wise = step(modelo_v2rmauto_ls, direction = c("both"), 
                            trace = F)


summary(modelo_v2rmauto_wise)
## 
## Call:
## lm(formula = V2RMAuto_std ~ Poblacion + ArCiclista + ArBajaVel + 
##     PMPeatones + PMCiclistas + PMTPublico + PMVMotor + Precipitacion + 
##     Temp, data = df_v2rmauto)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -36.676 -15.748  -2.734  10.980  32.301 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     62.263      4.693  13.266 2.55e-09 ***
## Poblacion       12.570      7.448   1.688  0.11361    
## ArCiclista     -14.850      6.127  -2.424  0.02950 *  
## ArBajaVel      -14.393      6.456  -2.229  0.04268 *  
## PMPeatones     469.573    180.098   2.607  0.02068 *  
## PMCiclistas    111.695     43.346   2.577  0.02195 *  
## PMTPublico     381.361    151.355   2.520  0.02452 *  
## PMVMotor       530.756    205.143   2.587  0.02150 *  
## Precipitacion   -8.022      6.936  -1.156  0.26684    
## Temp            19.378      6.163   3.144  0.00717 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 22.99 on 14 degrees of freedom
## Multiple R-squared:  0.7919, Adjusted R-squared:  0.6582 
## F-statistic: 5.921 on 9 and 14 DF,  p-value: 0.001701
barplot(log(vif(modelo_v2rmauto_wise))[vif(modelo_v2rmauto_wise) > 5],
          main = "log VIF DE LAS COVARIABLES", col = "blue")

abline(h = log(5), lwd = 2, col = "red")


METODO DE REGULARIZACION LASSO TARGET V2RMAUTO


xx = model.matrix(V2RMAuto_std ~., data = df_v2rmauto)[, -1]

modelo_v2rmauto_lassocv <- cv.glmnet(xx, y = df_v2rmauto$V2RMAuto_std, 
                                         alpha = 1,
                                         nfolds = 5)
modelo_v2rmauto_lassocv_min <- glmnet(xx, y = df_v2rmauto$V2RMAuto_std, 
                                         alpha = 1,
                                         lambda = modelo_v2rmauto_lassocv$lambda.min)

modelo_v2rmauto_lassocv_1sd <- glmnet(xx, y = df_v2rmauto$V2RMAuto_std, 
                                         alpha = 1,
                                         lambda = modelo_v2rmauto_lassocv$lambda.1se)

cbind(coef(modelo_v2rmauto_lassocv_min),
      coef(modelo_v2rmauto_lassocv_1sd))
## 12 x 2 sparse Matrix of class "dgCMatrix"
##                       s0       s0
## (Intercept)   62.2628935 62.26289
## Poblacion      0.6046567  .      
## DenPob         .          .      
## ArCiclista    -7.9051384  .      
## ArBajaVel      .          .      
## PMPeatones     .          .      
## PMCiclistas    .          .      
## PMTPublico     .          .      
## PMVMotor       .          .      
## Precipitacion -2.7152213  .      
## Temp          19.6160028 10.04084
## PBI            .          .

Vamos a adoptar el criterio min en este caso

modelo_v2rmauto_ls_2=lm(V2RMAuto_std ~ Poblacion + ArCiclista +
                        Precipitacion + Temp, data = df_v2rmauto)

summary(modelo_v2rmauto_ls_2)
## 
## Call:
## lm(formula = V2RMAuto_std ~ Poblacion + ArCiclista + Precipitacion + 
##     Temp, data = df_v2rmauto)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -57.942  -9.842  -1.198  11.807  54.438 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     62.263      5.354  11.630 4.39e-10 ***
## Poblacion        2.223      6.448   0.345  0.73411    
## ArCiclista     -14.727      5.897  -2.498  0.02185 *  
## Precipitacion   -7.668      7.201  -1.065  0.30032    
## Temp            22.959      6.397   3.589  0.00196 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 26.23 on 19 degrees of freedom
## Multiple R-squared:  0.6326, Adjusted R-squared:  0.5552 
## F-statistic: 8.178 on 4 and 19 DF,  p-value: 0.0005183
modelo_v2rmauto_wise_2=step(modelo_v2rmauto_ls_2, direction = c("both"), 
                            trace = F)

summary(modelo_v2rmauto_wise_2)
## 
## Call:
## lm(formula = V2RMAuto_std ~ ArCiclista + Precipitacion + Temp, 
##     data = df_v2rmauto)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -57.636 -11.680  -0.996  11.353  51.952 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     62.263      5.234  11.895 1.59e-10 ***
## ArCiclista     -15.379      5.461  -2.816  0.01067 *  
## Precipitacion   -8.731      6.362  -1.372  0.18515    
## Temp            23.043      6.250   3.687  0.00146 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 25.64 on 20 degrees of freedom
## Multiple R-squared:  0.6303, Adjusted R-squared:  0.5748 
## F-statistic: 11.37 on 3 and 20 DF,  p-value: 0.0001438
modelo_MM_v2rmauto= lmrobdetMM(V2RMAuto_std ~ Poblacion + ArCiclista + Precipitacion + 
    Temp, data = df_v2rmauto, control = control)
plot(modelo_MM_v2rmauto$residuals, modelo_MM_v2rmauto$rweights)

Observamos dos ciudades que son outliers

sort(modelo_MM_v2rmauto$rweights)[1:2]
## 6 9 
## 0 0
modelo_v2rmauto_wise_3=step(lm(V2RMAuto_std ~ ., data = df_v2rmauto[-c(6, 9),]), 
                            direction = c("both"), 
                            trace = F)

summary(modelo_v2rmauto_wise_3)
## 
## Call:
## lm(formula = V2RMAuto_std ~ Poblacion + PMVMotor + Temp + PBI, 
##     data = df_v2rmauto[-c(6, 9), ])
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -20.094  -6.085  -1.078   7.216  32.995 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   54.166      2.853  18.985 7.00e-13 ***
## Poblacion     18.051      3.257   5.543 3.58e-05 ***
## PMVMotor      10.612      3.108   3.415  0.00330 ** 
## Temp          11.825      3.242   3.647  0.00199 ** 
## PBI            3.958      2.918   1.356  0.19275    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 13.21 on 17 degrees of freedom
## Multiple R-squared:  0.8061, Adjusted R-squared:  0.7605 
## F-statistic: 17.67 on 4 and 17 DF,  p-value: 6.9e-06
par(mfrow=c(2,2))
plot(modelo_v2rmauto_wise_3)

No se observa estructura, ni residuos mayores a modulos de 3.



TARGET AUTOSM


modelo_autosm_ls=lm(AutoSM_std ~., data = df_autosm)

summary(modelo_autosm_ls)
## 
## Call:
## lm(formula = AutoSM_std ~ ., data = df_autosm)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -13.7315  -4.7622   0.5719   5.7051  16.7060 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    17.6636     2.1624   8.169 3.03e-06 ***
## Poblacion      -0.5314     3.5538  -0.150   0.8836    
## DenPob         -3.7350     3.9084  -0.956   0.3581    
## ArCiclista     -3.0996     2.8946  -1.071   0.3053    
## ArBajaVel      -5.9694     3.3017  -1.808   0.0957 .  
## PMPeatones     49.8306    88.3383   0.564   0.5831    
## PMCiclistas    12.6611    21.0795   0.601   0.5593    
## PMTPublico     47.7785    73.9904   0.646   0.5306    
## PMVMotor       61.0973   100.2855   0.609   0.5537    
## Precipitacion  -2.4117     3.2566  -0.741   0.4732    
## Temp            3.0742     3.3159   0.927   0.3721    
## PBI             2.1197     2.6883   0.788   0.4457    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 10.59 on 12 degrees of freedom
## Multiple R-squared:  0.663,  Adjusted R-squared:  0.3541 
## F-statistic: 2.146 on 11 and 12 DF,  p-value: 0.1026
barplot(log(vif(modelo_autosm_ls))[vif(modelo_autosm_ls) > 5],
          main = "log VIF DE LAS COVARIABLES", col = "blue")

abline(h = log(5), lwd = 2, col = "red")


METODO DE REGULARIZACION LASSO TARGET AUTOSM


xx = model.matrix(AutoSM_std ~., data = df_autosm)[, -1]

modelo_autosm_lassocv <- cv.glmnet(xx, y = df_autosm$AutoSM_std, 
                                         alpha = 1,
                                         nfolds = 5,
                                         maxit=1000000)
modelo_autosm_lassocv_min <- glmnet(xx, y = df_autosm$AutoSM_std, 
                                         alpha = 1,
                                         lambda = modelo_autosm_lassocv$lambda.min)

modelo_autosm_lassocv_1sd <- glmnet(xx, y = df_autosm$AutoSM_std, 
                                         alpha = 1,
                                         lambda = modelo_autosm_lassocv$lambda.1se)

cbind(coef(modelo_autosm_lassocv_min),
      coef(modelo_autosm_lassocv_1sd))
## 12 x 2 sparse Matrix of class "dgCMatrix"
##                      s0         s0
## (Intercept)   17.663610 17.6636102
## Poblacion      .         .        
## DenPob         .         .        
## ArCiclista    -1.763633 -0.5199541
## ArBajaVel     -4.210453 -2.9179391
## PMPeatones    -2.254634 -1.2171255
## PMCiclistas    .         .        
## PMTPublico     .         .        
## PMVMotor       .         .        
## Precipitacion  .         .        
## Temp           .         .        
## PBI            .         .
modelo_autosm_ls2 = lm(AutoSM_std ~ ArCiclista + ArBajaVel + PMPeatones, 
                       data = df_autosm)

summary(modelo_autosm_ls2)
## 
## Call:
## lm(formula = AutoSM_std ~ ArCiclista + ArBajaVel + PMPeatones, 
##     data = df_autosm)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -11.927  -6.676  -1.188   6.273  17.445 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   17.664      1.919   9.203 1.25e-08 ***
## ArCiclista    -3.427      2.160  -1.587   0.1282    
## ArBajaVel     -5.940      2.127  -2.793   0.0112 *  
## PMPeatones    -3.643      2.234  -1.631   0.1186    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.402 on 20 degrees of freedom
## Multiple R-squared:  0.5575, Adjusted R-squared:  0.4912 
## F-statistic: 8.401 on 3 and 20 DF,  p-value: 0.0008222
par(mfrow=c(2,2))
plot(modelo_autosm_ls2)

A partir de los graficos de los residuos, decidimos intentar alguna transformacion a ver si mejora nuestro modelo.

df_autosm$AutoSM_std = df_autosm$AutoSM_std + 0.01

boxcox(modelo_autosm_ls2)

df_autosm$AutoSM_std = sqrt(df_autosm$AutoSM_std)

modelo_autosm_ls3 = lm(AutoSM_std ~ ArCiclista + ArBajaVel + PMPeatones, 
                       data = df_autosm)

summary(modelo_autosm_ls3)
## 
## Call:
## lm(formula = AutoSM_std ~ ArCiclista + ArBajaVel + PMPeatones, 
##     data = df_autosm)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.86367 -0.74229  0.06008  1.03209  1.52277 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   3.8928     0.2300  16.922 2.55e-13 ***
## ArCiclista   -0.4568     0.2589  -1.765  0.09290 .  
## ArBajaVel    -0.7946     0.2549  -3.117  0.00543 ** 
## PMPeatones   -0.3729     0.2677  -1.393  0.17896    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.127 on 20 degrees of freedom
## Multiple R-squared:  0.5799, Adjusted R-squared:  0.5169 
## F-statistic: 9.203 on 3 and 20 DF,  p-value: 0.0004977
par(mfrow=c(2,2))
plot(modelo_autosm_ls3)

modelo_MM_autosm = lmrobdetMM(AutoSM_std ~ ArCiclista + ArBajaVel + PMPeatones, 
    data = df_autosm, control = control)
plot(modelo_MM_autosm$residuals, modelo_MM_autosm$rweights)

No vemos outliers a considerar



TARGET AUTOAUTO


modelo_autoauto_ls = lm(AutoAuto_std ~., data = df_autoauto)

summary(modelo_autoauto_ls)
## 
## Call:
## lm(formula = AutoAuto_std ~ ., data = df_autoauto)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.99348 -0.37553  0.00617  0.24130  2.51788 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)
## (Intercept)   -1.137e-15  1.958e-01   0.000    1.000
## Poblacion     -1.467e-01  3.218e-01  -0.456    0.657
## DenPob         1.595e-01  3.539e-01   0.451    0.660
## ArCiclista    -1.533e-01  2.621e-01  -0.585    0.570
## ArBajaVel     -2.961e-01  2.990e-01  -0.990    0.342
## PMPeatones     1.215e+00  7.999e+00   0.152    0.882
## PMCiclistas    2.535e-01  1.909e+00   0.133    0.897
## PMTPublico     1.536e+00  6.700e+00   0.229    0.823
## PMVMotor       1.965e+00  9.081e+00   0.216    0.832
## Precipitacion -3.464e-01  2.949e-01  -1.175    0.263
## Temp          -2.556e-01  3.003e-01  -0.851    0.411
## PBI           -7.630e-02  2.434e-01  -0.313    0.759
## 
## Residual standard error: 0.9593 on 12 degrees of freedom
## Multiple R-squared:  0.5199, Adjusted R-squared:  0.07983 
## F-statistic: 1.181 on 11 and 12 DF,  p-value: 0.3878

Realizamos un estudio de la colinealidad de las covariables del modelo, ya que tan pocos betas significativos nos da la idea de un problema por ese lado.

barplot(log(vif(modelo_autoauto_ls))[vif(modelo_autoauto_ls) > 5],
          main = "log VIF DE LAS COVARIABLES", col = "blue")

abline(h = log(5), lwd = 2, col = "red")


APLICAMOS METODO LASSO DE REGULARIZACION TARGET AUTOAUTO


xx = model.matrix(AutoAuto_std ~., data = df_autoauto)[, -1]

modelo_autoauto_lassocv <- cv.glmnet(xx, y = df_autoauto$AutoAuto_std, 
                                         alpha = 1,
                                         nfolds = 5)
modelo_autoauto_lassocv_min <- glmnet(xx, y = df_autoauto$AutoAuto_std, 
                                         alpha = 1,
                                         lambda = modelo_autoauto_lassocv$lambda.min)

modelo_autoauto_lassocv_1sd <- glmnet(xx, y = df_autoauto$AutoAuto_std, 
                                         alpha = 1,
                                         lambda = modelo_autoauto_lassocv$lambda.1se)

cbind(coef(modelo_autoauto_lassocv_min),
      coef(modelo_autoauto_lassocv_1sd))
## 12 x 2 sparse Matrix of class "dgCMatrix"
##                          s0            s0
## (Intercept)   -1.698226e-16 -1.457168e-16
## Poblacion      .             0.000000e+00
## DenPob         .             .           
## ArCiclista     .             .           
## ArBajaVel      .             .           
## PMPeatones    -2.895011e-01  .           
## PMCiclistas    .             .           
## PMTPublico     .             .           
## PMVMotor       .             .           
## Precipitacion  .             .           
## Temp           .             .           
## PBI            .             .

Si bien nuestro modelo lasso nos devuelve PMCiclistas y PMPeatones, vamos a optar por quedarnos solo con PMPeatones, ya que por construccion sabemos que ambas covariables tienen problemas de multicolinealidad elevada.

modelo_autoauto_ls_2 = lm(AutoAuto_std ~ ArCiclista + ArBajaVel + PMPeatones, data = df_autoauto)

summary(modelo_autoauto_ls_2)
## 
## Call:
## lm(formula = AutoAuto_std ~ ArCiclista + ArBajaVel + PMPeatones, 
##     data = df_autoauto)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.0828 -0.3846 -0.2230  0.2706  2.8322 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)  
## (Intercept) -1.919e-16  1.650e-01   0.000   1.0000  
## ArCiclista  -1.956e-01  1.857e-01  -1.053   0.3047  
## ArBajaVel   -1.712e-01  1.828e-01  -0.936   0.3602  
## PMPeatones  -4.647e-01  1.921e-01  -2.420   0.0252 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8084 on 20 degrees of freedom
## Multiple R-squared:  0.4318, Adjusted R-squared:  0.3465 
## F-statistic: 5.066 on 3 and 20 DF,  p-value: 0.009019
summary(step(modelo_autoauto_ls_2, direction = c("both"), trace = F))
## 
## Call:
## lm(formula = AutoAuto_std ~ PMPeatones, data = df_autoauto)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.0520 -0.5155 -0.2052  0.1780  2.9979 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)   
## (Intercept) -2.394e-16  1.662e-01   0.000  1.00000   
## PMPeatones  -6.049e-01  1.698e-01  -3.563  0.00174 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8142 on 22 degrees of freedom
## Multiple R-squared:  0.3659, Adjusted R-squared:  0.3371 
## F-statistic: 12.69 on 1 and 22 DF,  p-value: 0.00174
modelo_MM_autoauto = lmrobdetMM(AutoAuto_std ~ ArCiclista + ArBajaVel + 
                               PMPeatones, data = df_autoauto, control = control)
plot(modelo_MM_autoauto$residuals, modelo_MM_autoauto$rweights)

Observamos 2 ciudades con valores de residuos muy elevados

sort(modelo_MM_autoauto$rweights)[1:2]
##         22         13 
## 0.00000000 0.01822396
modelo_autoauto_ls_3=lm(AutoAuto_std ~ ArCiclista + ArBajaVel + PMPeatones, 
    data = df_autoauto[-c(13, 22),])

summary(modelo_autoauto_ls_3)
## 
## Call:
## lm(formula = AutoAuto_std ~ ArCiclista + ArBajaVel + PMPeatones, 
##     data = df_autoauto[-c(13, 22), ])
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.81049 -0.17920  0.05436  0.21974  0.60267 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept) -0.19813    0.07612  -2.603  0.01800 * 
## ArCiclista  -0.15063    0.08161  -1.846  0.08145 . 
## ArBajaVel   -0.07311    0.08097  -0.903  0.37848   
## PMPeatones  -0.27977    0.08756  -3.195  0.00502 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3547 on 18 degrees of freedom
## Multiple R-squared:  0.5857, Adjusted R-squared:  0.5167 
## F-statistic: 8.483 on 3 and 18 DF,  p-value: 0.001003
par(mfrow=c(2,2))
plot(modelo_autoauto_ls_3)

Observamos una buena apariencia en los residuos

cor.test(abs(modelo_autoauto_ls_3$residuals), 
         modelo_autoauto_ls_3$fitted.values, method = "spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  abs(modelo_autoauto_ls_3$residuals) and modelo_autoauto_ls_3$fitted.values
## S = 1532, p-value = 0.5478
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##      rho 
## 0.134952

Comprobamos que no hay Heterocedasticidad

vif(modelo_autoauto_ls_3)
## ArCiclista  ArBajaVel PMPeatones 
##   1.167712   1.128402   1.224928

Eliminamos el problema de multicolinealidad entre las covariables



Vamos ahora a analizar todas las clases de accidentes viales, utilizando la variable accidentes_viales, que es una combinacion lineal de las anteriores

TARGET ACCIDENTES VIALES


modelo_accidentes_viales_ls = lm(accidentes_viales_std ~. ,df_accidentesviales)

summary(modelo_accidentes_viales_ls)
## 
## Call:
## lm(formula = accidentes_viales_std ~ ., data = df_accidentesviales)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.18745 -0.44059  0.06232  0.43035  1.38822 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)  
## (Intercept)   -6.477e-15  1.582e-01   0.000   1.0000  
## Poblacion      2.280e-01  2.600e-01   0.877   0.3977  
## DenPob        -1.974e-02  2.859e-01  -0.069   0.9461  
## ArCiclista    -3.354e-01  2.118e-01  -1.584   0.1393  
## ArBajaVel     -4.172e-01  2.416e-01  -1.727   0.1097  
## PMPeatones     1.228e+01  6.463e+00   1.900   0.0817 .
## PMCiclistas    2.958e+00  1.542e+00   1.918   0.0792 .
## PMTPublico     1.064e+01  5.413e+00   1.965   0.0730 .
## PMVMotor       1.445e+01  7.337e+00   1.970   0.0724 .
## Precipitacion -2.692e-01  2.383e-01  -1.130   0.2807  
## Temp           1.019e-02  2.426e-01   0.042   0.9672  
## PBI           -4.025e-02  1.967e-01  -0.205   0.8413  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.775 on 12 degrees of freedom
## Multiple R-squared:  0.6866, Adjusted R-squared:  0.3993 
## F-statistic:  2.39 on 11 and 12 DF,  p-value: 0.07493
barplot(log(vif(modelo_accidentes_viales_ls))[vif(modelo_accidentes_viales_ls) > 5],
          main = "log VIF DE LAS COVARIABLES", col = "blue")

abline(h = log(5), lwd = 2, col = "red")


APLICAMOS METODO LASSO DE REGULARIZACION TARGET ACCIDENTES_VIALES


xx = model.matrix(accidentes_viales_std ~., data = df_accidentesviales)[, -1]

modelo_accidentesviales_lassocv <- cv.glmnet(xx, y = df_accidentesviales$accidentes_viales_std, 
                                             alpha = 1,
                                             nfolds = 5)
modelo_accidentesviales_lassocv_min <- glmnet(xx, y = df_accidentesviales$accidentes_viales_std, 
                                         alpha = 1,
                                         lambda = modelo_accidentesviales_lassocv$lambda.min)

modelo_accidentesviales_lassocv_1sd <- glmnet(xx, y = df_accidentesviales$accidentes_viales_std, 
                                         alpha = 1,
                                         lambda = modelo_accidentesviales_lassocv$lambda.1se)

cbind(coef(modelo_accidentesviales_lassocv_min),
      coef(modelo_accidentesviales_lassocv_1sd))
## 12 x 2 sparse Matrix of class "dgCMatrix"
##                          s0            s0
## (Intercept)   -1.901963e-16 -1.656208e-16
## Poblacion      1.906163e-02  .           
## DenPob         .             .           
## ArCiclista    -3.022070e-01 -1.261378e-01
## ArBajaVel     -2.195146e-02  .           
## PMPeatones    -2.311678e-01 -5.059179e-02
## PMCiclistas    .             .           
## PMTPublico     .             .           
## PMVMotor       .             .           
## Precipitacion  .             .           
## Temp           .             .           
## PBI            .             .
modelo_accidentes_viales_ls2 = lm(accidentes_viales_std ~ Poblacion + ArCiclista 
                                                        + ArBajaVel + PMPeatones, data = df_accidentesviales)

summary(modelo_accidentes_viales_ls2)
## 
## Call:
## lm(formula = accidentes_viales_std ~ Poblacion + ArCiclista + 
##     ArBajaVel + PMPeatones, data = df_accidentesviales)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.31256 -0.50546 -0.03908  0.40300  1.37683 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)  
## (Intercept)  1.846e-17  1.505e-01   0.000   1.0000  
## Poblacion    2.466e-01  1.631e-01   1.512   0.1471  
## ArCiclista  -3.563e-01  1.762e-01  -2.022   0.0575 .
## ArBajaVel   -2.046e-01  1.715e-01  -1.193   0.2474  
## PMPeatones  -3.516e-01  1.756e-01  -2.002   0.0598 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7371 on 19 degrees of freedom
## Multiple R-squared:  0.5511, Adjusted R-squared:  0.4566 
## F-statistic: 5.832 on 4 and 19 DF,  p-value: 0.00309

Observamos que los coeficientes estimados de nuestra regresion tienen un pvalue, cercano al 5%, deberemos tener en cuenta esto.

vif(modelo_accidentes_viales_ls2)
##  Poblacion ArCiclista  ArBajaVel PMPeatones 
##   1.126095   1.314217   1.244303   1.305666

No se observa problemas de multicolinealidad entre las covariables.

par(mfrow=c(2,2))
plot(modelo_accidentes_viales_ls2)

Observamos que la observacion 20, tiene un muy alto leverage y un residuo de 2, es una ciudad que se debe estudiar con cuidado. Fuera de eso, no se ve estructura en los residuos y parece haber homocedasticidad y distribucion normal de los residuos. En el QQplot, volvemos a ver resaltada la observacion 22.

modelo_MM_accidentesviales = lmrobdetMM(accidentes_viales_std ~ Poblacion + ArCiclista + 
    ArBajaVel + PMPeatones, data = df_accidentesviales, control = control)

plot(modelo_MM_accidentesviales$residuals, modelo_MM_accidentesviales$rweights)

No se encuentran outliers. No nos encontramos que se le haya asignado un rweight de 0 a la observacion 22, con lo que optamos por no considerarla outlier y la conservaremos.

summary(modelo_MM_accidentesviales)
## 
## Call:
## lmrobdetMM(formula = accidentes_viales_std ~ Poblacion + ArCiclista + ArBajaVel + 
##     PMPeatones, data = df_accidentesviales, control = control)
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -1.334816 -0.473743  0.008481  0.402913  1.446983 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.03804    0.17184  -0.221    0.827
## Poblacion    0.21387    0.19004   1.125    0.274
## ArCiclista  -0.30853    0.19437  -1.587    0.129
## ArBajaVel   -0.22743    0.19878  -1.144    0.267
## PMPeatones  -0.32415    0.19935  -1.626    0.120
## 
## Robust residual standard error: 0.8036 
## Multiple R-squared:  0.5982, Adjusted R-squared:  0.5136 
## Convergence in 37 IRWLS iterations

Si bien el R2adj de nuestro modelo robusto es superior al de OLS anteriormente mencionado (modelo_accidentes_viales_ls2), vamos a preferir el anterior para tareas de inferencia, debido a que sus coeficientes son mas significativos que este caso.